Response functions of cold neutron matter: density, spin and current fluctuations 
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We study the response of a single-component pair-correlated baryonic Fermi-liquid to density, 
spin, and their current perturbations. A complete set of response functions is derived in the low- 
temperature regime both within an effective theory based on a small momentum transfer expansion 
and within a numerical scheme valid for arbitrary momentum transfers. A comparison of these 
two approaches validates the perturbative approximation within the domain of its convergence. We 
derive the spectral functions of collective excitations associated with the density, density-current, 
spin, and spin-current perturbations. The dispersion relations of density and spin fluctuations are 
derived and it is shown that the density fluctuations lead to exciton-like undamped bound states, 
whereas the spin excitations correspond to diffusive modes above the pair-breaking threshold. The 
contribution of the collective pair-breaking modes to the specific heat of neutron matter at subnuclear 
densities is computed and is shown to be comparable to that of the degenerate electron gas at not 
too low temperatures. 

PACS numbers: 97.60. Jd,26.60.+c,21. 65. +f,13.15.+g 



I. INTRODUCTION 

The interiors of neutron stars become supcrnuid 
shortly after their formation (for reviews of the physics 
of superfluidity in neutron stars see Refs. [IMEl)- In the 
inner crust of a neutron star the neutrons pair in the 1 So 
channel with the density-dependent gap parameter in the 
range A < 1 MeV [J-Q. The neutron S-w&ve superflu- 
idity persists up to the densities of order of the nuclear 
saturation density no. Neutron P-wave superfluidity is 
expected at larger densities [6 8]. Protons, which are less 
abundant, form an S-wave pair condensate from densi- 
ties ~ no/2, where they de-confine from crustal nuclei, 
up to densities n no, i.e., the deep interiors of the 
star. For not too large isospin asymmetries, the D-wave 
condensation of neutron-proton pairs may set in at high 
densities as well @. 

The low-energy dynamics of baryonic matter in com- 
pact stars can be described, microscopically, in terms of 
a set of response functions to perturbations having dif- 
ferent symmetries. A frequently encountered example is 
the radiation and transport of neutrinos, in which case 
one is interested in vector and axial-vector perturbing 
operators. Response functions also contain the complete 
information on the spectrum of the low-lying excitations 
(i.e., density waves, spin waves, etc.) and, therefore, 
they permit to evaluate the contribution of the collective 
excitations to thermodynamics and transport of matter. 

Near equilibrium the response functions of nuclear sys- 
tems are characterized by length scales that are large 
compared to the inverse Fermi wave vector, or cquiva- 
lently, energies that are small compared to the Fermi en- 
ergy. In the unpaired limit, the Landau theory of normal 
Fermi liquids provides a suitable framework for the evalu- 
ation of response functions in compact stars [Iol - [T6| . The 
many-body problem of the evaluation of response func- 
tions entails a number of challenges. One is the identifi- 
cation of the relevant set of diagrams, when perturbation 
theory fails. For many systems the response functions are 



computed from a resummation of an infinite number of 
finite temperature ring diagrams [l7| . While this scheme 
accounts for the (vertex renormalized) single particle- 
hole excitations, it does not include multi-pair contribu- 
tions to the response functions. Such contributions are 
important for the evaluation of the magnetic suscepti- 
bility of degenerate nuclear matter [U>tl4| . The second 
challenge is the inclusion of the non-central forces, which 
arise in the nuclear systems due to the tensor forces. In 
fact, these give rise to the coupling of states with more 
than one quasiparticle-quasihole pair, thereby changing 
the static susceptibility and the magnetic moments in 
the nuclear Fermi liquid [IoHT3 |. As a consequence, in 
nuclear matter the relationship between the Landau pa- 
rameters and the magnetic susceptibility is considerably 
more complicated than for systems with purely central 
forces. Sum-rule arguments can be used to place a lower 
bound on the contribution to the static susceptibility 
coming from transitions to multipair states [ToL llIT | . Fur- 
thermore, it was shown that the rates of processes in- 
volving transitions to two quasiparticle-quasihole states 
may be calculated in terms of the collision integral in 
the Landau transport equation for quasiparticles [r2,ll3j]. 
The multi-loop processes induced by the tensor forces are 
of paramount importance in the astrophysics of neutron 
stars, since the bremsstrahlung processes on weak neutral 
currents are among the leading processes contributing to 
the neutrino luminosity of these stars [l8l - l2i| . 

The focus of this paper is the derivation of the re- 
sponse functions associated with perturbations of den- 
sity, density current, spin, and spin current in a single- 
component Fermi liquid. It extends our earlier study of 
density response [22j to new types of perturbations as 
well as revises some of the perturbative results contained 
therein. The energy scale characterizing the dynamical 
processes in neutron stars are of the order of temperature 
T < 1 MeV. The high densities in compact stars render 
the Fermi energies of fermions in the range cf ~ 10—100 
MeV. Consequently, one needs the response functions in 
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the limit T/ep <C 1. As mentioned above the pairing gaps 
could be of the order of 1 MeV, i.e., they substantially 
influence the dynamics of the systems for temperatures 
T < T c , the critical temperature of the supernuid phase 
transition. 

This study is based on the method of the Green's 
functions for supcrfluid systems at non-zero tempera- 
tures and aims at the resummation of an infinite se- 
ries of particle-hole ladder diagrams in neutron matter. 
This re-summation scheme respects the gauge invariance, 
sum rules, and baryon number conservation. The appro- 
priate technique was first developed by Abrikosov and 
Gor'kov in the electrodynamics of superconductors [23| 
(see also Ref. [HI)- I n this theory the response of the 
superconductors to external probes is expressed in the 
language of propagators at non-zero temperature and 
density with contact interactions that do not distinguish 
among the particle-hole and particle-particle channels. 
It is equivalent to the theories initially advanced by Bo- 
golyubov (25[, Anderson (2(| and others, which are based 
on the equations of motion for second-quantized opera- 
tors. Subsequently, Larkin, Migdal, and Leggett (27l.[28j 
generalized the Landau Fermi-liquid theory to supercon- 
ductors and superfluids, thus extending the Abrikosov- 
Gor'kov approach to strongly interacting regime. This 
last method implements the wave-function renormaliza- 
tion of the quasiparticle spectrum, higher order harmon- 
ics in the interaction channels, and postulates particle- 
hole (ph) and particle-particle (pp) interactions with dif- 
ferent strength and/or sign. 

The response functions of baryonic matter were studied 
in the unpaired, but degenerate regime in the context of 
neutrino emission from compact stars (see, e.g., Ref. [29[ 
and references therein). The work on these functions 
in the same context, but for supcrfluid baryonic matter 
started more recently [30M37j . 

Quite generally, the response functions to density, spin 
and their current perturbations can be related to the ap- 
propriate response functions of baryons to the operators 
of the electrowcak theory. To see the mapping explicitly 
consider the weak interaction Lagrangian, which at low 
energies is given by 

where Gf is the Fermi constant and the vector and axial- 
vector currents are defined as 

J£ = c v ^ n Y^n 2tc v ip i N (l,v F )il> N , (2) 
J'X = CA*JV7 M 75*JV - ca^Pn {(tv f ,ct) ip N , (3) 

and = V'7 A i(l — 75)^ is the lepton current. Here v F 
is the Fermi velocity of baryons, <x is the vector of Pauli- 
matrices. Here and below the Greek indices run over 0,1, 
2, 3, and label the temporal and three spatial coordinates; 
the spatial coordinates are also labeled by Latin indices 
and run through 1, 2, 3. Equations © and ^ approxi- 
mate the baryonic vector and axial-vector weak currents 



by their dominant contributions in the non-relativistic 
limit by keeping the large components of the baryonic 
Dirac spinors. The bare vertices of interest are thus given 
by the expression in-between the baryon fields ip N and 
ip N in Eqs. © and ©: 



r?" = (r°,r°) = (i,v F ), (4) 
= (r s ,r s ) = (<Tv F ,cT). (5) 



It is now clear that there is a one-to-one correspon- 
dence between weak interaction vertices in the non- 
relativistic limit and vertices associated with the density 
and density-current (index D), as well as the spin-current 
and spin-density perturbations (index S). 

One complication that is always present in compact 
stars is the fact that the matter is multi-component in 
the crusts and the core of the star. A supernuid fea- 
tures Goldstone bosons associated with the breaking of 
the baryon U(l) number in a supernuid [38l [391 ] . Fur- 
thermore, the existence of the lattice of nuclei (and non- 
spherical nuclear phases) in the crust adds the lattice 
phonons to the set of the collective modes that propa- 
gate in the star's crust [4(| EH- The various modes are 
coupled [39|, |42T|44| . In the cores of neutron stars there 
are at least three fluids — the neutron and proton Fermi 
liquids, which are both expected to be in the supcrfluid 
state, and an ultra- relativistic gas of electrons [20]. The 
density modes associated with the superconducting pro- 
ton component in the homogeneous matter of the outer 
core of neutron stars were computed in Refs. j45l - |47| . It 
is clear that our treatment of a single-component super- 
fluid nuclear Fermi liquid does not account for coupling 
among various components. A more complete treatment 
must take into account the multi-component nature of 
matter. 

This paper is organized as follows. The remainder of 
the Introduction provides prerequisite information. In 
Sec. [XT] the baryon propagators and self-energies are in- 
troduced within a finite-temperature imaginary-time the- 
ory. Vertex functions corresponding to density and spin 
perturbations are discussed in Sec. IIII1 Section flV] is de- 
voted to the density and spin response functions, with 
two subsections discussing perturbativc expansions of re- 
sponse functions as well as their exact numerical evalua- 
tion. In Sec. [V] the spectral functions and collective den- 
sity and spin excitations are discussed. We evaluate the 
specific heat contribution arising from these excitations 
in Sec. lVIl Our conclusions are collected in Sec. I VIII The 
details of computations are relegated to Appendices E] 
and [B] and a comparison to other methods is presented 
in Appendix [Cj We use the natural units fi = c = 1 and 
assume that the Boltzmann constant ks = 1, with the 
exception of Sec. IVII 
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A. Prerequisites 

In this study we explore the temperature domain well 
below the critical temperature of supcrfluid phase transi- 
tions; typically T/T c < 0.5, where T c is the critical tem- 
perature of a supcrfluid transition. This is the case in 
the dominant majority of observable neutron stars. We 
further consider densities where the 1 So-wave pairing is 
dominant among neutrons and protons. This assumption 
confines our study to the densities at and below the nu- 
clear saturation density. In the presumed temperature 
and density domain it is safe to treat the nucleons as 
non-relativistic particles, i.e., the Fermi velocity of the 
particles is small compared to the velocity of light in a 
vacuum, vp <C 1 in natural units. This enables us to 
use the non-relativistic dispersion law for the particles in 
the normal state and non-relativistic limits of the Dirac 
matrices appearing in the bare vertices. Furthermore, be- 
cause we work in the extreme low-temperature limit, we 
shall restrict the length of the momenta of the particles to 
their Fermi wave- vector, i.e., we write p = m*vpn, where 
to* is the effective mass of a quasiparticlc and n = p/\p\. 

One of the purposes of this work is to compare the re- 
sponse functions obtained from perturbative approaches 
and direct numerical computation. The perturbative 
treatment is based on a low-momentum transfer expan- 
sion, where the expansion parameter is either generic 
and reflects the characteristic properties of the system 
or is dictated by certain kincmatical conditions valid in 
the domain of interest. Examples of small parameters 
arc q/kp or qvpjw, where hp is the Fermi wave vec- 
tor, and u) and q are the energy and the magnitude of 
the momentum transfer. While the first parameter is 
generic for thermal processes (i.e., processes in which 
the energy-momentum transfer is of the order of temper- 
ature) the second is small only in the kincmatical domain 
of time-like processes (e.g., neutrino radiation). In the 
second case the momentum transfer is thermal, therefore 
q <C hp, which establishes one suitable expansion param- 
eter. We note that for on-shell perturbations with linear 
spectrum, as, for example, neutrinos (u = q in natural 
units) the smallncss of the two expansion parameters re- 
duces to the condition vp <C 1, which is the same as the 
non-relativistic expansion. In the case of the numerical 
computation there are in principle no constraints on the 
values of the momentum transfer and the Fermi wave 
vector. However, since our intention is to compare the 
perturbative and exact numerical results, we will restrict 
ourselves to the range of values of the parameters defined 
by the perturbative treatment. 

B. Unpaired and pair-correlated particle spectra 

As we work in the non-relativistic limit the spectrum 
in the normal state is given by 

Cp = h * (6) 



where \i is the chemical potential. The spectrum in the 
pair-correlated case is 

e p = y/% + A2, (7) 

where we assume that the gap function is momentum 
independent, which is the case for contact pairing inter- 
actions. We will need frequently the perturbed spectra 
of particles, which arc defined in the unpaired 

where in the second expression the small recoil term 
q 2 /8m* has been dropped. In the paired case the quasi- 
particlc spectrum is 

e± = y/& + A 2 ^ yfej ± Z P qv, (9) 
to leading order in |q|. 

II. BARYON PROPAGATORS AND 
SELF-ENERGIES 

In a normal Fermi liquid the propagator is defined as 
G n ,^(p,t-t') = -8 aal {T T i> pa (T)^ pa ,{r')), (10) 

where r is the imaginary time, a is the spin projection, 
and T T is the time-ordering operator. The Dyson equa- 
tion for the normal propagator is given by 

Gn = Gq + GoSGat, (11) 

where the index refers to the free-particle propagator 
and £ is the self-energy. A superfluid is described by the 
following propagators 

G aa ,{p,T-T') = -5 CTCT '(T T ^ CT (r)^(r')), (12) 
F aa ,(p,T- T ') = (T T V>- pi (r)VH(0), (13) 
F+,( P ,t-t') = (T T ^ pt (r)^_ pi (r')), (14) 
G- a ,{p,T-r') = -5^{T T il pa {T)^ pa ,{T')). 

(15) 

These propagators obey Nambu-Gorkov equations and 
are given by 

G = G q + G q £g + G AF + =G N +G N AF+, (16) 
F+ = G„£-F+ + G„A + G = G N A+G, (17) 

F = GoSF + GoAG" = GatAG", (18) 
G~ = Gq +G^t-G- + GoA + F = G N + G N A+F, 

(19) 

where G(p) and Go(p) are the full and free normal prop- 
agators, F(p) and F + (p) are the anomalous propagators, 
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£(p) and S (p) are the normal self-energies for particles 
and holes, and A(p) and A + (p) are the anomalous self- 
energies. The propagators and self-energies are 2x2- 
matrices in the spin space. (We suppress the isospin 
space variables as we consider only single-component en- 
sembles with fixed isospin.) The normal (particle and 
hole) propagators and self-energies are diagonal in spin 
space, 

G(p) = G(p)l 2 = G"(-p)i a = G-(-p), (20) 
£(p) = S(p)i 2 = E-(-p)l 2 = S-(-p), (21) 

while the anomalous ones are antisymmetric in spin space 
and therefore are proportional to ia 2 , 

F(p) = F(p)ia 2 , F+(p)=F+(p)ia 2 (22) 

A(p) - A(p)ia 2 , A + (p)ia 2 =A + (p)ia 2 , (23) 

where a 2 stands for the second Pauli matrix. For real 
pairing gaps A(p) = A + (p) and F(p) = F + (p). 

The propagators can be written as the sum of a pole 
and a regular part by expanding the self-energy in the 
vicinity of the Fermi surface. Neglecting the (small) off- 
shell contributions, we shall keep the pole part of the 
propagators and set the wave function rcnormalization 
Zip) -1 = 1 — cL£(k-0U=£ = 1- The real-time solution 
of the Nambu-Gorkov equations in momentum space are 
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with 



Pa - £ P + ii] Po + £ P + if] 

G(p) = G-(-p), 

E(p) = E-(-p), 

F(p) = F + {p)=F{p), 

A( P ) = A+(p) = A, 

and the Bogolyubov amplitudes defined as 

1 (. t v " 



V2 



1 



(24) 



, (25) 



(26) 
(27) 
(28) 
(29) 



(30) 



(31) 



The finite temperature Matsubara Green's functions are 
obtained via a replacement of the time-component of the 
four-momentum in Eqs. (|24]l and (|25|) by a complex fre- 
quency 



which assumes discrete values p n = (2n + 1)itT, where n 
is an integer. 



III. VERTEX FUNCTIONS 

The equations for vertex functions involve loops which 
are constructed from the convolutions of a product of 
two propagators. One possible kinematics for such prod- 
ucts is the symmetrical one, which assigns to an arbi- 
trary imaginary-time propagator X the arguments X+ = 
[ip n + iu> m ,p+ |) and X_ = (i Pn ,p - f), i.e., the ex- 
ternal momentum is split symmetrically among the par- 
ticle and the hole (but the energy transfer is not). The 
remainder of this work will use this kinematics. We now 
turn to the calculation of the effective (or dressed) ver- 
tices, which take into account the modifications due to 
the strong interactions in the medium. The driving inter- 
action in the particle-particle and particle-hole channel 
will be parametrized as [27j 



(34) 
(35) 



where V D and V s are the interaction strengths in the 
density and spin channels, the subscripts or superscripts 
pp and ph refer to the particle-particle and particle-hole 
channels, respectively. 

Since the particle momenta are restricted to the Fermi 
surfaces, the amplitudes will depend only on the angle 
formed by the momenta of the particles. Therefore, as in 
the ordinary Fermi-liquid theory, they can be expanded 
in spherical harmonics with respect to this angle. The 
coefficients in this expansion arc the Landau parameters. 
We will retain the leading-order Landau parameter only, 
since the higher-order Landau parameters are numeri- 
cally insignificant. We will use below their values for 
bulk neutron matter as computed in Ref. [33j . 

In analogy with the random phase approximation for 
unpaired ensembles the calculation of full vertices re- 
quires a summation of an infinitely long chain of irre- 
ducible particle-hole ring diagrams. One possible way 
to derive these equations is to compute the variations of 
the Nambu-Gor'kov equations in an external field [27j . 
Another method to set up the integral equations for the 
vertices is to construct them directly from Feynman di- 
agrammatic rules. In any case, since a single-component 
supcrfluid ensemble is fully described by four different 
propagators, one finds that there are four topologically 
different vertices, which are determined by four coupled 
integral equations. The analytical form of these equa- 
tions for scalar vertices is 



G(ip n ,p) 
F(iPn,p) 



V 



iPn - £ 



>P,, 
1 



1 



(32) 
(33) 



r 



D/S 



\D/S 



+FT° /S G + GT^'°F + FVY /0 F 



d 4 p 



[2n)H 



V 



D/S 



ph 



D/S t 



Gff /S G 



;D/S } 



(36) 
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d 4 p 



\D/S 



(27T) 
+ ( 



Gr? /S F + Fts /S F}, (37) 



+Fff /s G + G-r 4 

(27r) 4 i 



M3/S 



D/S 



F + FT 



D/S-, 



H.D/S- _ 







r>/s 



F ,(38) 



(39) 



where subscripts -D and S* refer to the density and spin, 
f is the bare vertex for holes. Identical equations 

can be written for vector vertices. In the following we 
approximate the particle-hole and particle-particle inter- 
action amplitudes by the leading-order Landau parame- 
ters v p h and v pp . The last of these is determined by the 
gap equation as follows 



1 



vv r _ 



A l-2/( £p ) 
Kv 2e„ ' 



(40) 



where v = m*kF/2ir 2 is the density of states on the Fermi 
surface and A is the cut-off which regularizes the ultra- 
violet divergence of the integral. 

The solutions of the vertex equations (|36[) to (|3U)) arc 
described in appendix ??. We find for bare scalar vertex 
1 



1 



rf(w,g)=rf(a; jQ ) = 



C(u,,q)-vP h Q+(uj,q) 

V + (u,q) 
C(u,q)-vP h Q+(u,q) 



, (41) 
(42) 



for the bare vector vertex TP = v 



rf /4 (w,9) 



±n„ 



VpiQ-^q) 



C( U ,q)-v^Q+(uj,q) 



(43) 



2/3 



± 



V-(u,q) 



+ 



C(co,q)-vP h Q+(u;,q) 



n q v F 

CM)' 



(44) 



for the bare scalar spin-current vertex Tq = crv 

v^ h Q+(cj,q)n q 



rf/ 4 (w,g) 



n„ ± 



C(w,q) - v° b Q-(u},q 



v F , 



(45) 



■ 2/3 



(u>,q) = ± 



V+(u,q) 



< h P-(^g)Q + (a;,g) 
C(u),q)-vS h Q-(uj,q) 



crn q vp 

cRgp 



and, finally, for the bare spin vertex Tq = a 

<rC(w,q) 



C{u,q)-v^Q-{u,q) 

(TV-(uj,q) 
' C(u,q)-v$ h Q-(uj,q) 



(46) 

(47) 
(48) 



Tf{u,q) = -Ti(u,q): 
rf(w,g) = -rf(w,g) 



The functions on the right-hand side of Eqs. (|4T[) to 
are defined in Appendix [X] The full vertex entering the 
density response is seen to coincide with the one derived 
in Refs. [22, HH, [Hj] . The remainder vertices are in agree- 
ment with the ones obtained in Ref. 13611 . 



IV. RESPONSE FUNCTIONS 

We start with a general expression for a response func- 
tion in terms of a current-current correlation function 



n ,u, = 1 

2 



d 4 p 
(2^ 



Tr 



(49) 



where Jq and J v are the bare and dressed currents. The 
polarization tensor consists of four different contributions 
(we drop here the subscripts D/S) 



2 
1 

+ 2 
1 

+ 2 
1 

+ 2 



d 4 p 
(2ir) A i 

d 4 p 
(2tt)H 

d A p 
(2tt) 4 £ 

d 4 p 
(2tt) 4 £ 



Tr 



Tr 



Tr 



Tr 



q 



'f^( P +|)^(p-| 



(50) 



The trace should be carried out in the spin space. We 
can now compute the response functions by substituting 
the bare and effective vertices corresponding to the de- 
sired type of perturbation. For the density response the 



vertices are Tq D and T® D and we find 



Q + (uJ,q) 



C(u,q)-v? h Q+(u;,q) 



(51) 



Furthermore, the density-current response is given by 
(summation over repeated indices is assumed) 



A-(u,q) 



B{LO,q)f>-(uj,q) 
C{u,q) 



v» h Q + (u,,q)Q-(u,,q) 
~C(cj,q) 2 -v? h C(uj,q)Q+(LJ,q) 



(52) 
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the spin-current response is given by 



ti$(u>,q) = <A + (u,q)- 



+ - 



v S ph Q-{u,q)Q+{u,q) 



C(u,q) 3 -V^C{u,q)Q-(u,q) ^ 
and finally, the spin-density response is 

3fi-(w,«) 



C(u,q) 

(53) 



(54) 



The functions appearing on the right-hand side of 
Eqs. ([5"2"jl to are defined in Appendix [X] For den- 
sity perturbations the off-diagonal elements of the polar- 
ization tensor with mixed temporal and spatial indices 
are given by (below for the sake of brevity we drop the 
arguments of the loops) 

dn \ A+C - BV+ 

- — s-^r) n v v F, (55) 



and 



v°(A+C-BV+)q 



ri> } v F , (56) 



while for spin-perturbations they are given by 

C J 

-nl }v F , (57) 



m?( w ,q) = J^U+ni-B (Sin; 

v s ph (A-C-BV- 



o,-, , fdQ [a~C-BV- 1 . . 

*s^ = ]^{ c- v ? h Q- ) n ^- (58) 

Each of the polarization tensors can be decomposed into 
transverse and longitudinal parts with respect to the di- 
rection of the momentum transfer q according to 



U L {u,q) = n 00 (u,q) 
1 



(59) 



Il T (u>,q) = -[8» -n»njjn«(w,g). (60) 

Performing the decomposition of the vector polarization 
tensor we obtain for the longitudinal projection 



ny,L(w,g) 



Q + (w,q) 



C(u,q)-v° h Q+(oj,q) 



(61) 



and for transverse projection 
n ViT (o;, Q ) = ^ J^}.A-(u,q)-A-(u,q)( 



n{n{f 



(62) 



The longitudinal and transverse components of the axial- 
vector polarization read 



U AiL (u,q) = lA + (u,q)- 



B(u>,q)T> + (u,q) 



C(u,q) 



+ 



ru,r(w,g) 



yS h Q-(u,q)Q+(u;,q) 
C(co, q y-v^ h C(u J ,q)Q-(u J ,q) 

Q-(u,q) 

C{w,q)-v^Q-{w,q) ' 



v% , (63) 



(64) 



These results, which arc valid for arbitrary orientations 
of the external vectors fields, can be further simplified by 
a suitable choice of the coordinate system. 



A. Perturbative results 

We now expand the loop functions with respect to the 
small parameter y = q/kp and keep contributions up to 
fourth order in this parameter. The thermal function Q 
depends on y and therefore, we can write 



fe=0 



G2k x Zk y 2k = Go + g 2 x"y 2 + x*y* + 0(y b ). 



(65) 

The expansions of the loop functions contain only even 
functions of the parameter y, since possible odd terms 
will disappear after angle integration; thus, e.g., for the 
,4-loop we obtain 



A = A +A 2 y 2 +A±y 4 + 0(y e ), 



(66) 



and similarly for the other three. In practice, we ex- 
pand the pre-factors in Eqs. (|A25[) to (|A28[) as well as 
the function Q in the power series in parameter y and 
subsequently combine them. This leads us to the follow- 
ing explicit expressions: 



-in 



4//V , U^X 4 4 

— a-» + . a y 



G(v,qv,q) 



Go+ {^ Go + l G2 ) y2 

16 M 4 A/1 2 i \ 



(67) 



= - / ^ 



dn 

4~ 



4:fi 2 X 2 2 Wfj, 4 x 4 4 



V r 2 



16/x 



. "2/ 

A ,,2 



g(v,qv,q) 



L V r 1 4 
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where we have dropped terms 0(y 5 ) and higher. Note 
that the term Q(q,v,qv) is purely real, i.e., does not 
contribute to the imaginary parts of the loops. For fixed 
momentum transfer it is constant and yields numerically 
negligible contribution. For the remaining loops we ob- 
tain 



v- x C = 
+ 

l v- = 



w r 


w 2 


^ a 


(_ M 2 


4A 2 Go ' 


I 3A 2 


f- 


Qi + 


\ 5A 2 


2 20A< 


2A Go + 


^ r 2 

^G 2 y 4 



10A 



G±y\ (69) 



Go 



12A 2 



10A 



G 4 y\ 



I' 



dCl 

4-7T 



-^Goy + G 3 y 



= 0. 



(70) 
(71) 

(72) 



One can now readily identify the coefficients of the ex- 
pansion (|66[) and its counterparts for the remaining loops. 
In full analogy, an expansion of the polarization tensors 
is given as 



ll D/S 



LL D/S,0 



nr - 4 



1 D/5,2J/" + llS/S.4y" + °(j/ )■ ( 73 ) 

The coefficients of the density response function are 
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ttOO 
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(74) 
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where 

Qo 



Co-v° h Q 



C 2 - v ph (A+C + A+C 2 - B 2 V+ - B Q V+ 

C 4 - v ph (A+C + A+C 2 + A+C 4 
B 4 V+ -B 2 V+ -B V+). 



(76) 

(77) 
I, 

(78) 
(79) 



However, it turns out that Qo = and 

„4+C 2 - B 2 CoV+ + BoC 2 V+ - BoC V+ = 0. (80) 

Consequently, the expansion coefficients of the polariza- 
tion tensor arc 



n oo 

n D0 
ttOO 
ii D2 

ttOO 
ii D4 



(81) 
(82) 



(83) 



i.e., the density response function obtains a non-zero con- 
tribution at order y 4 . The coefficients for the current 
response are 



0, 
0, 

i [A+Co + A+C 2 + A^C 4 

<--0 \ 

B 4 V+ - B 2 V+ - BoV+] , 



TT" - 



it" = < .At 

L1 D2 ~ \ -^2 



(84) 



li D4 



Aa 



Vph 



BxD7 v° h AT Co {A+ Co - B{D+ 



BxCqVT - B X C 2 VT + B 3 C VT 

c 2 



F I 



(85) 



(A3 Cl + 2A\CpC 2 - Ar<#) (AjCp - B{D+ 
C 2 



where 



AT L4+C 2 + A+Co - B X V+ - B 3 V^ 



2C C 2 - v° h [A+Cl + A+C C 2 - B0C0V+ 
B C 2 V+ - B 2 CoV+ 



■F- 

(86) 



(87) 



In the case of the current response the first non-zero term 
arises at the order y 2 and the fourth order term is sub- 
leading. 

Note that the vector current polarization tensor must 
vanish at the zeroth order as required by the /-sum 
rule 17 



lim / dej u> lmIl D (q , uj) = 0. 
J 



(88) 



This is a direct consequence of the conservation of the 
baryon number. For the spin-current response we find 



TT UU - A + n z 

11 SO — -^0 V Fi 



(89) 



n 5 ° 2 = {At + 1& {A+ATCoAiATCo )\v 2 F , (90) 



n°° 4 = \At + vUAiA-Co + A+A~Co + A+A7C 2 
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-A^B,V+ 



A^B 3 V+ 



- A^B X V+ 
AtBxA^Co 



J F- 



(91) 



In this case the leading order contribution, given by 
Eq. (|89l) . is of order y° in the y-expansion, which means 
that the vertex corrections introduce sub-leading order 
corrections and the single-loop result is a good approxi- 
mation to the full polarization tensor. Finally, the spin- 
density response is given by 



TT3J 

TT-" 

il S4 



0. 

3{A± 



') V 



(92) 
(93) 
(94) 



The leading order contribution now arises at order y 2 . 

If we restrict ourselves only to the leading order con- 
tributions in each channel, then these contain only the 
leading order term in the expansion of the thermal func- 
tion C/Qj i- e. , at order y°. The explicit expressions are 



ttOO _ 
11 D — 


6 V 4r 
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~ 45w« 


(95) 
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(96) 


ttOO _ 
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IT" - 


-^y Go=- 


9 9 

qv F r 


(98) 



The last equalities in these expressions make it clear that 
the expansion, which was initially carried out with re- 
spect to the parameter q/kp maps onto the expansion 
in Vp. It is seen that the vector current polarization 
tensors arc of order 0(v F ) while the axial vector polar- 
ization tensors are of order 0(v F ). The pcrturbative re- 
sults (p?Sl and (|M|) are in good agreement with the ones 
derived recently in the context of vector neutrino emis- 
sion (31T 13(1 |37| . Similarly, the perturbative expressions 
in the spin channel ([§?]) and ([§8"]) are in agreement with 
the original results derived in the context of the axial 
vector neutrino emission [30, |48h50I1 . 



B. Numerical results for response functions 

Figures Q] and [5] show the dependence of the real and 
imaginary parts of the density and spin response func- 
tions, respectively, on the transferred energy for fixed 
three-momentum transfer. The zero temperature gap is 
fixed at A(0) = 1 MeV and T c = A(0)/1.76. The low- 
est order Landau parameter is set v® h = —0.5 for den- 
sity perturbations and v ph — 0.5 for spin perturbations 
(these correspond to the values computed in Ref. [H|). 
The frequency and momentum transfer are normalized to 
the threshold frequency 2A(T). The response function in 
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Figure 1: (Color online) Numerical (solid lines) and pertur- 
bative (dashed lines) results for the imaginary (heavy, blue 
line) and the real (light, cyan line) parts of the (a) density 
response function ilj^ and (b) the density-current response 
function 11^ normalized to the density of states v in units 
fin -2 . The energy transfer uj is in units of 2A(T). The tem- 
perature is fixed at 0.5 T c with pairing gap A = 1.0 MeV. 
The ratio of momentum transfer and Fermi momentum is 
kept fixed at q/fci? = 0.01 and the Fermi momentum is set to 
= 1.0 fin - , which translates to the density n = 0.221no. 



the negative energy range can be obtained from the re- 
lations ReIV D / s {-Lo) = Rcn c / S (w) and ImlF 2 '< ' s '(-w) = 
—ImIl D > s (uj). The numerical method of computing the 
response functions exactly is described in appendix ??. 

Our comparison of the perturbative analytical results 
with the exact numerical ones shows that (i) for the 
density response the higher-order corrections shift the 
imaginary part to higher frequencies, i.e., for a fixed fre- 
quency the imaginary part is larger; the real parts are 
correspondingly larger as well, (ii) For the density cur- 
rent response the perturbative and exact results match 
to a high accuracy; (iii) for the spin-current response 
both results match again to a high accuracy; (iv) for the 
spin-density response small deviations arc observed close 
to the threshold; the imaginary part is again shifted to 
higher frequencies. Note that in each case the imaginary 




parts arc identically zero below the threshold for pair 
breaking process 2A(T). 

The density response function can be compared to the 
one derived in a previous paper [22j . As shown above, 
the first non-vanishing contribution arises from the term 
n™ 4 and not from IT™ as in Ref. 0, where II™ ^ 0. 
Consequently, the numerical values of the real and imagi- 
nary parts are roughly by an order of magnitude smaller. 
Nevertheless, the dependence of the real and imaginary 
parts of the polarization tensor on the frequency shows 
essentially the same behavior. The difference between the 



present results and that of Ref. [22j can be understood 
as follows. We note that the general form of the density 
response function in [22j , Eq. (18) and the definitions of 
the elementary loops, Eqs. (19) to (22), are the same. 
The difference arises at the level of the loops A, B, and 
C given by Eq. (25) to (27) of [Hf. The most general 
form of the first loop, upon substitution of Bogolyubov 
amplitudes in Eqs. (19) and (21) of Ref. [22[ is given by 



A{q) 



d 3 p 
(2^)3 



[( e + e ')( ee ' - ££' + A 2 ) + u(£'e - £e') 



^ (99) 

with the short-hand notations £ = S, p - q /2, C = £ P +q/2 




1.03 



Figure 3: (Color online) The spectral function of (a) the den- 
sity fluctuation, (b) current fluctuations and (c) spin-density 
fluctuations as a function of the energy transfer u and of 
momentum transfer in units of 2A(T) at T = 0.5T C and 
1:f = 1 fin - , which corresponds to density 0.211 no, where 
no = 0.16 fm -3 is the nuclear saturation density. 
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and e = y/? + A 2 , e' = V^') 2 + A 2 and 



given by the solution u>(q) of the equation 



1 

2tl' 



ui- 



(e + e'f 



(100) 



l-^ h ReP(o; ) q) = 0. 



(105) 



We see that Eq. (25) of Ref. [22| does not contain the 
w(£'e — £e') which vanishes manifestly in the limit q — > 0, 
but is finite if q ^ 0. For the remaining B and C loops 
we obtain 



B(q) 
C(q) 



2A 



d 3 p 



(2tt) 3 
d 3 p fl-2/(e) 




we' + (e' + e)£']Sf (101) 

(e + e')(«' + + A 2 ) 

(102) 



and we see that the terms (e' + e)£' and — U)(^e' + £'e) are 
missing in Eqs. (26) and (27) of Ref. [22|. Both terms 
that were dropped vanish in the limit q — > 0, because 
they arc odd in £, while after changing the integration 
measure according to Eq. (|A24[) , we obtain integrals over 
symmetrical in £ limits. Thus, we conclude that, the 
discrepancy between the present treatment and that of 
Ref. [22J originates from incomplete expressions in Eqs. 
(25)- (27) of the latter work. 



V. SPECTRAL FUNCTIONS AND 
COLLECTIVE MODES 

The knowledge of the response functions allows us to 
construct an effective theory of excitations in the nu- 
clear medium. Their full (interacting) propagator is com- 
pletely determined by their spectral function, which in 
each channel is defined via the imaginary part of the po- 
larization as R(uj,q) = — 2lmH(u, q). For example, in 
the density channel, using Eq. (|5Tj) , one finds 



R(u,q) = - 



2(^)- 2 ImP( W ,q) 



(103) 



where P(u>,q) = Q + (ui , q) / C (uj , q) . Similar relations hold 
for other excitation channels (i.e., current-density, spin- 
current and spin-density). Above the threshold u> > 
2A(T) non-zero imaginary part implies that the collec- 
tive excitations have finite life-time, i.e., are not per- 
fect quasiparticles. Nevertheless, in the limit where the 
imaginary part is small one can approximate the spectral 
function as 

R(u>, q) = 2TrZ(q)6((v° h )- 1 - ReP(uj, q)) + R TCg (uj, q), 

(104) 

where R IOS (uj,q) is the regular (i.e. smooth) part of the 
spectral functions and Z(q) is the wave-function renor- 
malization. The dispersion relation of the excitations is 



Figure [3] shows the dependence of the spectral func- 
tions for density [Fig.[3[a)], current [Fig. [3jb)] , and spin- 
density [Fig. [22c)] fluctuations on the energy and momen- 
tum transfer. The spectral functions have a Breit-Wigner 
form, therefore the location of their maxima is controlled 
by the real parts of the response functions, whereas their 
widths by the imaginary parts. In the case of density 
fluctuations the imaginary component of the polariza- 
tion tensor has a power-law (oc g 4 ) behavior for fixed en- 
ergy transfer, as is explicit from the analytical form (|95|) . 
At fixed momentum transfer the spectral function has a 
threshold due to the proportionality Qo oc $( w — 2A). 
At low-momentum transfers the main contribution to 
the spectral function comes from the vicinity of the pair 
breaking threshold (uj ~ 2A); for large momentum trans- 
fers, modes away from the energy threshold become im- 
portant. The qualitative features seen in the spectral 
function of the density response are seen also for the cur- 
rent response; some quantitive differences arise because 
now, for fixed energy transfer, the imaginary part scales 
as oc q 2 , c.f. Eq. (|96[) . Consequently, the low-momentum 
contributions are only weakly suppressed and the maxi- 
mum of the spectral function is numerically larger. The 
response functions associated with spin perturbations ap- 
pear at order vp, therefore their absolute scale is larger 
than that for the density and current-density responses, 
which scale as v F . It has the same functional dependence 
on the momentum and energy transfer as the density- 
current response [see Eq. ([98)) ] and differs only by the 
numerical pre-factor and the vp dependence. For small 
momentum transfers the main contribution to the spec- 
tral function comes from the region near the threshold. 
Note that the spin-current response, to leading order, 
is independent of the momentum transfer, therefore the 
two-dimensional form given in Fig. H](c) is sufficient. Its 
dependence on the frequency reflects the dependence of 
the function Q , c.f. Eq. f9"7). 

From the spectral functions we can extract the quasi- 
particlc spectra of the collective excitations. These can 
be defined by the poles of the spectral function when 
lmP(ui,q) = 0, i.e., by the condition (|105|) . We start by 
setting the parameters characterizing the superfluid state 
to their relevant scales and by studying the nature of the 
modes as a function of the particle-hole interaction v p i L . 
The stability of the normal Fermi-liquid state constrains 
Vph > ~ 1 (note that we work at leading order in the ex- 
pansion of Landau parameters in spherical harmonics). 
The numerical solutions of Eq. (|105[) for the density and 
spin excitations are shown in Fig.2]for — 1 < w p h < 2 and 
fixed k F = 1 fin" 1 , T/T c = 0.5 with A = 1 MeV. For 
positive values of the particle-hole interaction the modes 
appear in the domain w/2A > 1, where lmP(ui,q) ^ 0, 
i.e., they represent damped (diffusive) modes of oscilla- 
tions of density and spin-density, respectively, associated 
with the pair-breaking processes. For negative values of 
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Figure 4: (Color online) Dispersion relations of collective excitations for density (left panel) and spin density (right panel) 
perturbations for the values of the particle-hole interaction v p h shown in each panel for /cf = 1 fm _1 and T/T c = 0.5. The 
heavy lines (blue online) correspond to undamped exitonic modes, the light lines (cyan online) correspond to diffusive damped 
pair-breaking modes. 



the particle-hole interaction, the modes exist in the do- 
main w/2A < 1, where the pair-breaking part of the 
ImP(w, q) vanishes; therefore the modes represent un- 
damped oscillations of density and spin-density around 
their average values. These modes are "exitonic" as they 
correspond to bound pairs of particles and holes. The 

Table I: The values of the coefficients in the fit formula (|106p 
for several values of the particle-hole interaction in the density 
(upper part) and spin (lower part) channels. The remaining 
parameters are fixed to their characteristic values JtF = 1 
fm _ , A = 1 MeV, and m*/m = 1. The interactions are 
given in units of the density of states u(kp). 





a 


b 


c 


d 


-1 


0.817879 


88.4029 


- 11411.2 




-0.5 


0.849483 


53.7703 


- 5210.88 




0.5 


0.600277 


80.4847 


115.199 




1 


0.586293 


57.3376 


49.0417 




2 


0.584029 


116.109 


193.221 




-1 


1.02616 


-84.0295 


2216.43 


-7449470. 


-0.5 


1.02487 


-39.6023 


-568.605 


-844334. 


0.5 


0.885174 


83.0049 


- 922.6 


4886.17 


1 


0.928416 


154.6 


-2912.19 


26000.1 


2 


0.924387 


313.859 


-12031.8 


218234. 



spectra in each case can be accurately fitted by the poly- 
nomial of the form 

Q(q) =a + bq 2 +cq i + dq 6 , (106) 

where in the case of density perturbations accurate re- 
sults are obtained with only three parameters (d = 0). 
Here we defined dimensionless quantities Cj = w/2A and 
q = q/kp. The fitted values of the parameters for the 
results shown in Fig. |4] and are given in Table Q] 

Next we consider a specific microscopic calcula- 
tion (33d , which provides us with the density dependence 



Table II: The same as in Table [J but for the density- 
dependent Uph interactions taken from Ref. [33| . The entries 
are the Fermi wave vector, the values of the coefficients in the 
fit formula (|106(l . and the range of momentum transfers Aq 
in units of Uf- 

}if t'ph a b c kp Aq 



-0.45 0.941366 2.69439 -34.0444 (0.197; 0.29) 

0.41 0.518545 12.4875 8.8972 (0.223; 0.3) 

-0.43 0.896992 12.9105 -447.747 (0.111; 0.22) 

0.40 0.6465 37.092 -197.621 (0.127; 0.3) 

-0.41 0.901629 67.8547 -12952.2 (0.047; 0.098) 

0.40 1.17054 84.3504 -488.556 (0.054; 0.3) 

-0.36 0.894476 702.973 -1283650 (0.015; 0.031) 

0.39 1.15211 859.001 -49727.6 (0.018; 0.08) 



of the parameters of the neutron superfluid and the as- 
sociated values of the leading-order Landau parameters 
in the particle-hole channel. We solved Eq. (|105[) in the 
density and spin channels for each density and subse- 
quently fitted the spectra with the formula (|106[) . The 
results are displayed in Table [TIJ some of the character- 
istics of the superfluid are shown in Table IIIII The den- 
sity excitations exist below the pair-breaking threshold, 
i.e., represent undamped exitonic modes. Conversely, be- 
cause Vph changes the sign in the spin channel, the spin 
excitations represent diffusive modes with finite damping. 
Note that each of these modes exist within some finite 
interval Aq of momentum transfers. The lower bound 
arises because perturbations that are sufficiently large to 
excite a mode arise at some finite value of q. The upper 
bound in most cases is the consequence of the use of per- 
turbativc response functions, whose validity breaks down 
for large momentum transfers q/kp ~ 0.3; in some cases 
the upper bounds are associated with the disappearance 
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of the solutions from the search domain. 



VI. SPECIFIC HEAT 

The specific heat contribution arising from the col- 
lective modes in the neutron star crust has recently at- 
tracted recently attention in the context of non-spherical 
phases [Jlj]. These modes at not too low temperatures 
can dominate the specific heat provided by the degener- 
ate, ultra-relativistic electron gas. Below we shall exam- 
ine the contribution of the collective modes discussed in 
the previous section to the specific heat of a supcrfluid 
neutron star crust. 

The entropy of a collective bosonic mode is given by 



S = 2k B ]T [(1 + g q ) ln(l + g q ) + g q 



In 



9q\ 



(107) 



Table III: Density and wave- vector dependence of the specific 
heat of various components in neutron matter and the crust 
of a neutron star at T — 0.5T C . Tabulated are the net den- 
sity of matter n n [fm -3 ], the neutron wave-vector fcF[fm -1 ], 
the neutron effective mass in units of bare mass, the pair- 
ing gap A[MeV], the electron wave-vector fcf e [fm _1 ], assum- 
ing that the electron number density n e = 0.023n n , the spe- 
cific heats of electron gas , Bogolyubov- Anderson phonons 



„(BA) 



pair-breaking density fluctuations Cy , pair-breaking 
spin fluctuations Cy in units of 10 18 erg cm -3 K _1 . 



fin 


Uf 


m* 


A k Fe 


c (e) 


„(BA) 

c v 




c (<t) 


0.034 


1.00 


0.94 


3.09 0.29 


16.6 


26.419 


2.479 


0.291 


0.058 


1.20 


0.92 


2.44 3.42 


18.9 


7.611 


3.508 


0.079 


0.093 


1.40 


0.88 


1.41 3.99 


14.8 


1.003 


0.008 


0.004 


0.138 


1.60 


0.84 


0.57 0.45 


7.8 


0.045 


0.012 


0.000 



where g q = [exp(oj q /T) — 1] _1 is the Bose distribution 
function of collective excitations with the spectrum u q . 
The specific heat is then given by 



dS 
kBT df 



\ - dg g 
2^ t dT' 



(108) 



For a collective (acoustic) mode with linear spectrum lo = 
u\q\, where u is the sound velocity, Eq. (|108[) can be 
integrated (55| 



>> 



2ir 2 k%T 3 



(109) 



An acoustic mode, in a compact star setting, is associ- 
ated with the nuclear lattice in the crust, where phonons 
contribute to the specific heat below the melting temper- 
ature of the crust ~ 10 9 K. At low temperatures the su- 
pcrfluid supports the Bogolyubov-Anderson (BA) mode 
with the velocity 



.BA 



V F 

7T 



^(i + ^ h ) 1/2 - 



(110) 



where vf = http/m* is the (effective) Fermi velocity. The 
dispersion relation (|110j) does not contain temperature 
corrections. In the following we will ignore the damping 
of the BA mode and extrapolate the result (|110j) to higher 
temperatures. Apart from these two collective modes, 
the main contribution to the specific heat of matter is 
due to the electrons which, in a first approximation, can 
be treated as a uniform ultra-relativistic ideal Fermi gas. 
At low temperatures their specific heat is then given by 



u 2 T 



3(hc) 3 



(111) 



where fi e is the electron chemical potential. 

Table Mil compares the various contributions to the 
specific heat of matter at subnuclear densities. The tem- 
perature at each density corresponds to T = 0.5T C , with 



T c = A/1.76. The contribution of the BA mode, c£ A is 
computed from Eqs. (|109[) and (|110j) . the contribution of 
electrons from Eq. (|111|) assuming n e = 0.023n n , where 
n e and n n are the electron and neutron number densities. 
The contributions from density and spin pair-breaking 
contributions, Cy and Cy ^ are computed through the nu- 
merical integration of Eq. (|108[) with the collective mode 
spectrum given by Eq. (|106|) . The coefficients a, b and 
c in Eq. (|106|) for the density fluctuations and the spin 
fluctuations, as well as the integration limits in Eq. (|108[) 
are tabulated in Table |TT] (the coefficient d = in all 
cases). It is seen that the density fluctuations consider- 
ably contribute to the net specific heat of matter for lower 
densities (wave- vectors) , while the spin- fluctuations are 
negligible at T/T c = 0.5. The result of for the BA mode 
should be taken as suggestive, because we neglected the 
temperature correction to the dispersion relation and the 
possible damping of this mode. 

The temperature dependence of the specific heat due 
to the pair-breaking modes and the specific heat of elec- 
tron gas is shown in Fig.[5]for kp = 1 fm _1 . The electron 
specific heat is linear in temperature, whereas the specific 
heat of the pair-breaking fluctuations has a power law be- 
havior, which is close to the T 3 law characteristic for lin- 
ear in q spectra. The difference reflects the non-linearity 
of the spectrum (|106|) . We have assumed that the tem- 
perature dependence of the coefficients a, 6, and c can 
be neglected in a first approximation, i.e., the spectrum 
of collective excitations is assumed to be independent of 
temperature. This assumption is validated by the insen- 
sitivity of the maxima of the spectral functions to the 
temperature variations (see Fig. [3]) which were compared 
at T/T c = 0.5 and 0.9. 



VII. CONCLUSIONS 

In this work we studied the response functions of a sin- 
gle component pair-correlated baryonic matter to den- 
sity, spin and their current perturbations in the low- 
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Figure 5: (Color online) Dependence of the specific heats 
due to electrons (solid line) and pair-breaking density (dashed 
line) and spin fluctuations (dash-dotted line) on the reduced 
temperature T/T c for kp = 1 fin - . 



temperature regime. These results should be relevant 
for the description of both the dynamical and thermody- 
namical properties of baryonic matter at low densities, 
i.e., the densities where the baryons form an S- wave su- 
pcrfluid. It was observed that the expansions in the pa- 
rameters q/kF and vfi/u lead essentially to the same 
perturbative results, which in turn can be interpreted 
as an expansion in the parameter vf <C 1- We have ap- 
plied an exact numerical method to evaluate the response 
functions and to validate the perturbative approximation 
in the domain of its convergence. We further derived the 
dispersion relations of the collective excitations of density 
and spin-density perturbations. For positive values of 
the particlc-hole interactions these correspond to weakly 
damped diffusive excitations, whereas for negative values 
- to undamped excitonic modes. 

The spectral functions presented above can be modi- 
fied in a number of ways. As noted in the Introduction 
the multi-loop processes were found to be important for 
the neutrino emission and they could additionally con- 
tribute to the spectral functions in the kinematical do- 
main where two-particle-two-hole excitations are impor- 
tant. Furthermore, higher order Landau parameters, if 
included into driving interactions in the particle-particle 
and particle-hole channels may re quir e some renormal- 
ization of the spectra, see Refs. |5ll - l54| . 

The application of the formalism to compute the 
specific heat of the matter expected in neutron star 
crusts shows that the contribution of the collective pair- 
breaking excitations can be a significant part of the net 
specific heat of matter. For some density parameters and 
not too low temperatures the combined contribution from 
superfluid modes of neutron fluid can be larger than the 
specific heat stored in the degenerate electron gas. 
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Appendix A: Solving the equations for the vertices 

The bare vertices given by Eq. ((4]) are diagonal in spin 
space. Likewise, the particle and hole vertices are diag- 
onal in spin space, i.e., Ti = T1I2 and F4 = T4l2. The 
anomalous vertices are proportional to the second Pauli 
matrix, T2 = Yiioi and = T^ia2- The equation for 
the hole vertex T4 can be obtained from the equation for 
the particle vertex Ti by interchanging particle and hole 
lines. To account for this property one can introduce, fol- 
lowing Ref. [53] , an operator V to revert the direction of 
ingoing and outgoing momenta and to exchange the spin 
indices simultaneously, when acting on a vertex function. 
The explicit action of this operator is 



rr 



l,a/3, 



(Al) 



whereby T is the time-reversal operator, i.e., it is equal 
to +1 for vertices which are even under time reversal 
operation and —1 for vertices which are odd under this 
transformation. By considering the action of the operator 
V on the bare vertices one finds that scalar vertices, e.g., 



Tq = 1 and Tq = crv do not change their sign, while 
vector vertices like Tq = v or Tq = er gain an additional 
minus sign. Furthermore, we note that the equations for 
T2 and T3 are adjoint to each other. Formally, one can 
cast this property into the equation 



r 



-r 2 



(A2) 



Note that the full current vertices can depend on any 
external momentum involved in the problem, therefore 
they need to be decomposed in components along the 
vectors v and q. Thus, the most general Ansatz for the 
density current vertices is 



*■ q "'qi 



r?n„ 



(A3) 
(A4) 



where the subscripts on the unit vector n refer to the 
vector defining its direction. The coefficients T® q are 
normalized such that T D n v = vp, i.e., the coefficient T D 
is simply the modulus of the Fermi velocity. Similar to 
Eqs. (|A3j) and (|A4[) decompositions holds for spin-current 
vertices. 

The solution of the system (f3"6"|) to (131)1) is simplified if 
one takes into account the identities 123, uSf 



GG'(v,LO,q) 
G-F(v,u,q) 
GF(v,tJ,q) 
G-G-(v,u,q) 



G-G{v,LO,q), 
-FG~{v,uj,q) 
-FG(v,Lo,q), 
GG(v,-u),q), 



(A5) 
(A6) 
(A7) 
(A8) 
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where the products of the Green's functions refer to their 
convolutions defined as z 



oo °? 

X+X'_ = vT I d£, p X ( 



ip n +iu) m ,p+ - 



72 — — OO 



X X' (ipn,p- 



(A9) 

where X± <G {G±,G±,F±,F±}. The solution contains 
the following linear combinations of the convolutions (241 

A(V) = G+G- - F+F^V, (A10) 
B = G+F_ - F+G^, (All) 
C = G+GZ+F+F--V-', (A12) 

D(V) = -G+F_ - F+GZV. (A13) 

We take the matrix structure of vertices and propagators 
into account and use the relations (jAip and (|A2[) to cast 
the set of the four coupled integral equations into the 
following two equations for the new vertex functions 

r-D/Si \ -nD/S, x , [ d£l ^ r D/S , / \ 

T ' a (v,u,q) = T/ (t>,o;,g) + / — V ph ' (v,v) 
x U(P)r £, / s '(u / ,a;, Q ) + Bf D / s {v',u,q)} , (A14) 

x "(C + Vp)f D/S (^,^ 9 ) + J D(P)r^ s ( 1 ;',c, g )" . 

(A15) 

To write down the solutions of the integral equation we 
need the following angle averages of the loop functions 



where dn = sin 9 dd dtp. Furthermore, we need the angle 
averages of first moments of the loop functions with re- 
spect to the cosine of the angle enclosed by n v and n q , 
i.e., x = n q ■ n v , which we write as 



y 



dn~ 

4ir ' 



(A20) 



where 



Y = xvT d£ P x( 

n=—oc 
(iPn,P 



ip n + lLU m ,p+ - 



x X' 



— OO 

2 



(A21) 



We also define the auxiliary combination of the loops: 

fi ± (w ) g) = A ± (oj, q)C(uj, q) - B(u, gJX^w, g),(A22) 
Q^u, g) = A ± (oj, g)C(w, g) - B(u, q)V + (uj, g).(A23) 

The computations of the phase-space integrals in 
Eqs. (jATOj) to ([A13l) can be simplified [28| , because each 
loop can be written as a product of some thermal func- 
tion and a pre-factor that depends only on the quantities 
u! and q ■ vp. To carry out the phase-space integrations 
we first change the integration measure: 



d 3 p 
(2tt) 3 



dn 

4^ / d ^ 



(A24) 



where we used the fact that at low temperatures the lower 
integration limit —fi/T ~ — oo. For the sake of complete- 
ness we list the resulting expressions for the loops [23, [HI 



A d 



dn 

47T 

qv 



l±V 



G(v,u>,q) 



B 
C 
V ± 



uj — qv 

' dil uj + qv 



4tt 2A 



dCl> 



-in 

dn 

47T 



g(v,qv,q) - G(v,uj,q 
V 



(qvf 



4A 2 

uj + qv uj — qv 



(A25) 

(A26) 
(A27) 



4A 



4A 



Q(v,u,q), 

(A28) 



A(V) = 


f dn 


(A16) 




B = 


r dn 

1 47T 


(A17) 




C = 




(A18) 




V(P) = 




(A19) 





where the thermal function is given by 



e+~e- /(£_)-/(£+) 
e+e_ cj 2 — (e+ — £-) 2 + if] 

e+ + e- l-/(e_)-/(e + ) 



e_i_e_ or 



(A29) 



where /(ar) = {exp[(cc — n)/T] + 1} _1 is the fermionic 
distribution function. In the following we focus on the 
pair-breaking part of Eq. (|A29[) given by 



-t-oo 

g ph (v,cj,q) = -A 2 Jd^ 



(f+ + e_) 
e+e- 

— OO 

l-/(6_)-/( £+ ) 

w 2 - (e+ + e_) 2 + w?' 



(A30) 



which is the dominant part of the response in the low- 
temperature domain. 
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Appendix B: Thermal function 
1. Analytical result 

Here we determine the real and imaginary parts of the 
zeroth order coefficient in the expansion of the thermal 
function. The first step is to use the generalized Dirac 
identity 



dzf(z) 



(z-( + ir]) n+1 



dzf(z) 
(z - C)" +1 



dzf(z)—S(z-C) (Bl) 

and to decompose the complex function at hand into real 
and imaginary parts. The imaginary part can be inte- 
grated analytically using the partial integration in the 
formula 



fdzf(z)s^(z-o = (-ir^l 



VCe [a,0]. 
(B2) 

Once the imaginary part is calculated, the real part can 
be obtained via the Kramers-Kronig relation 



Re0(u;) = - 

7T 



+00 

1 f du' lm<j)(uj') 



(B3) 



provided that the imaginary part decays faster than 1/w 
for large w. The application of this procedure to the 
thermal function to leading order gives 



c7 p 



P b 



-A 2 



2tanh(^) 



4e 2 



i5) 



(B4) 



Next we compute the imaginary part and obtain 



Img o pb 



ttA 2 J dti p 

—00 
27rA 2 tanh 



2tanh „. , . 
S(uj 2 - 4e 2 ) 



,4T) 



\u>\ Vu 2 - 4A 2 V2 



^-A 



(B5) 



Note the threshold behavior enforced by the Heavyside 
function: Energy transfer is possible only for frequen- 
cies larger than the pair-breaking threshold 2A. Fur- 
thermore, the thermal function at this order is inde- 
pendent of the momentum transfer; this implies that 
the momentum-transfer dependence of the response func- 
tions is determined by the pre-factors of the loop func- 
tions [c.f. Eqs. (fAlB] to (|Al8|) ]. 



2. Numerical calculation of the thermal function 

In this section we focus on the numerical calculation of 
the angle average of the thermal function, which is given 



by 



d 3 p (e++e_) 



)3 £+£ _ 

(l-/(e_)-/(e + )) 
{u 2 - (e+ + e_) 2 + irj) 

27T +1 OO 

-2A 2 J— J— I 

-1 



,/e 2 -A 2 e + e - 



(l-/(6_)-/( e+ )) 

oj 2 - (e + + e_) 2 + if] 



(B6) 



The factor of 2 in the second relation arises from the fact 
that the integrand is an even function of £ p and the in- 
tegration can be restricted to the positive values of the 
argument. We have also replaced the integration over 
the unpaired spectrum by the integration over the paired 
spectrum by means of the relation = e p de p . The 

integral over the azimuthal angle is trivial, since the inte- 
grand is independent of <f>. Further, after using the Dirac 
identity we obtain 

+1 00 
dx f de p e p (e + + e_) 



Img* h (v,u;, q ) = 2nA 2 J^-J 



A V 4~ ^ t+e - 



x (1 - / (e+) - / (6_)) 5 (lo 2 (e+ + e_) 2 ) . (B7) 

One of the integrations can be carried out with the help 
of the S function. One finds 



TTA 2 



dtp £p 



E 



a y " p 

(e+ + e_)[l-/(e+)-/(e_)] 



e+e_ 



9(1- \x 1>2 \). (B8) 



where are the solutions of the equation cj — 

[e + (x) + £- (a;)] 2 = and the prime denotes a deriva- 
tive with respect to x; its explicit form is given else- 
where [531 ]. Once the imaginary part is computed, the 
real part follows from the Kramers-Kronig relation. This 
completes our numerical procedure for computing the re- 
sponse functions. Each of the loops (|A25|) to (|A28[) can 
be computed by multiplying the numerical result for the 
thermal function by the appropriate pre-factor. 



Appendix C: Comparison to Leggett's formalism 

Here we compare the response functions derived above 
with the results obtained in the Leggett formalism [28|, 
HH and establish the correspondence between the two. 
The full (effective) normal vertices in the Leggett's 
formalism [28j j are defined as symmetrical and anti- 
symmtcrical combinations of the particle r and hole T h 
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vertices 



r = -{t + t 



(CI) 



If r = ±t , i.e., the vertices have odd or even parity 
under transformations which convert particles into holes, 
then one of the linear combinations (|C1|) vanishes. The 
bare vertices are defined analogously 



s- 



r-i(e-«* 



(C2) 



The anomalous vertex is denoted by f. With these defi- 
nitions the integral equations for the full vertices are 

1-/^0"') a 



4tt 



2A 2 



V 

4tt pp 



V^ s (vv')^X(v')r+ 



0. 



(C3) 



dfi' ..d/s, ,, qv' . , . 

(„„ )-^-A(w )r 



r 



1 



4tt P h 



4^ V (™)^«(«)r + =£ 
4^V (™)^7^K 



(C4) 



1 



^V^ S (vv')( K (v')-2X(v') 



(C5) 



where 



«(«') 


1 

2 




G+G— 


+ F+F_, (C6) 


A(t/) 




-F-, 






(C7) 






lim 


GG h + 




(C8) 













where, as before, the wave- function renormalization is set 
to unity. Keeping only the lowest order term in the ex- 
pansion of the particle-particle interaction in Eq. (|C3[) . 
one finds 1 — V^ p Ao = 0, i.e., the first two terms in that 
equation mutually cancel. (Note the different sign con- 
vention for Vp° p in the main body of the paper.) 

We proceed now to solve these equations for the ver- 
tices in some cases of interest. For that purpose define 
the following integrals: 



(C9) 



T) = 


r dn' 

1 4^T 




)«(»'), 


(CIO) 


7 = 


r <m' 

1 ~4t7 




) cos 2 6X(v') , 


(Cll) 


P = 


r dn' 

1 TtT 


V ph > (vv 


) cos 2 6 n(v') , 


(C12) 


i' = 


r dn' 

1 4^T 


v ( vv 


cos8 X(v') , 


(C13) 


= 


r dn' 

1 4^T 


V ph ' (vv 


) cos' 1 6 X(v') . 


(C14) 



In the following, we keep (as in the main body of this pa- 
per) the leading order Landau parameter in the particle- 
hole interaction amplitude, i.e., (vv') — u p h- Be- 
cause k and A are even functions of cos 9, in this ap- 
proximation the functions <j> and ip vanish. Following 
Leggett (28[, we will use below the abbreviations 



UJ 



qV F 
A ' 



(C15) 



The longitudinal component in the vector channel, is 
obtained when £ + = 1 and £~ = 0. Equations (|C3j) to 
(IC5I) are then written as 



^ s 2 a — 7 2uip —2 sua 

v ph uip 1 - v ph 77 v ph s (j) 

\ -u p h su a Uph s (f> 1 - w p h (17 - 2a) 



T~ 
T+ 





I 



As stated above = tp = at leading order in the 
particle-hole interaction. The solution of this matrix 
equation is given by 



2A s 2 c 
ui (s 2 a — 7UI 
1 



Vph 1 



1 - «ph Ql 



and r =0, where 



2 cry 



(C16) 
(C17) 

(C18) 



(s 2 a — 7) ' 

in agreement with Eqs. (52), (53) and (54) of Ref. [UJ 
taken in the case of V^ h = 0. 

The longitudinal projection of the vector current po- 
larization tensor is given by the expression [cf. [28j ]. Eq. 
(23a)] 



n 



V,L 



dn 

47T 



UJ UJ 

-rXf + —At" + (k- 2A) t 
A qVp 



(C19) 

which after the substitution of the vertices becomes 



n 



V,L 



2aj 



(s 2 a — 7) 



1 



1 - u P h Ql 1 - « P h Ql 

(C20) 
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By matching this equation to our result given by Eq. ([61 
we find 



V+ Q+ 



(C21) 



The latter equality is straightforward to prove by noting 
that Eqs. (|C7[) and (|C6[) can be written in terms of the 
thermal function (|A29|) as 



A = ~Q(v,u),q), 
f {qvf 



uj 2 - {qvf 



G(v,qv,q) - G(v,u,q) 



(C22) 



(C23) 



We conclude that our result for the longitudinal vec- 
tor current response function agrees with those given in 
Refs. [1^, [3^. In particular, we have verified that the 
limiting cases of (i) v p h = 0, (ii) cj <C A, qvp <C A 
for non-zero T, and (iii) same as in (ii), but for T = 0, 
we recover the results of Ref. [35| by using the match- 
ing condition (|C21[) . However, the perturbative result 
for the imaginary part of the longitudinal vector current 
response function in Ref. [H| [second term in Eq. (82)] 
differs from our result, given by Eqs. (p?5j) and (|B5[) by a 
factor of 1/8. (Note that the author Ref. [35[ used a den- 
sity of state which is by a factor of 2 larger than ours). 
We have verified that one recovers our result by starting 
from the exact expression (|C18|) for Ql and expanding 
the functions a, r\ and 7 in small s _1 . In Ref. [35j the 
exact expression is first approximated by Qj, ~ -i] + 2s~ 2r y 
after which the expansions for rj and 7 are substituted. 
The first step is the source of the discrepancy; we have 
verified that if the expansions of the functions a, rj and 
7 are directly substituted in the exact expression (|C18[) 
for Ql, then one recovers our result, which is also in 
agreement with the one quoted earlier by the authors of 

Refs. mm. 



In the case of the current response the bare vertices 
are given by £ = v± and £' 1 = — v±, i.e., 



r = o, r = »x 



(C24) 



where v± denotes the transverse to the momentum 
transfer projection of the quasiparticle velocity. If the 
momentum transfer is along the z axis, then v±_ = 
«i?(sin 9 cos <p, sin 9 sin <p, cos 9) . Keeping only the leading 
order Landau parameter in the particle-hole channel one 
finds 



0, r+«=0, T-^=vf> 



(C25) 



The anomalous vertex vanishes identically (the contri- 
butions to the vertex in the direction of the momentum 
transfer are neglected here). After substituting the ver- 
tices into the expression for the transverse part of the 
polarization tensor (Eq. (37) in Ref. [Hj]) we find 



n 



V,T 



VF 

2 



"k(1 



■(v-0), (C26) 



I h 

which coincides with Eq. (86) of Ref. [35| when ir h = 0. 
The transverse vector response function is given accord- 
ing to Eq. (|62p above. After substituting the explicit ex- 
pression for the loop function A~ from Eq. (|A25[) we use 
the relation (|C23j) to recover Eq. (|C26[) . i.e., the trans- 
verse vector polarization tensors are the same in both ap- 
proaches. However, the perturbative expansions of these 
transverse polarization tensors differ, by a factor 0{\), 
c.f. Eq. (88) in Ref. ^ and Eqs. JMJ and flE} above. 

The longitudinal and transverse axial-vector current 
polarization tensors of Ref. [35| can be matched to our 
results as in the case of vector current response functions, 
therefore we do not repeat the arguments above. 
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